function f=KCalculation(i,j,KaiA)
%Calculation of Kij

Khalf=0.43;
kxy0=[0 0 0 0.21;
      0.12 0 0.31 0;
      0 0 0 0.11;
      0 0 0 0];      %matrix(4,4)
kxyA=[0 0.332923 0 0.0798462;
      0.123 0 -0.319385 0;
      0 0.505692 0 -0.133077;
      0.479077 0 0.0532308 0];   %matrix(4,4)
  
f=kxy0(i,j)+kxyA(i,j)*KaiA/(Khalf+KaiA);  
return;